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SUMMARY 

A new and simple method of finite-element grid improvement is presented. 
The objective is to improve the accuracy of the analysis. The procedure is 
based on a minimization of the trace of the stiffness matrix. For a broad 
class of problems this minimization is seen to be equivalent to minimizing the 
potential energy. The method is illustrated with the classical tapered bar 
problem examined earlier by Prager and by Masur. Identical results are 
obtained. 


INTRODUCTION 

In a general context, the finite-element method is an approximate proce- 
dure for solving differential equations. The accuracy of the method depends 
on (1) the number of elements, (2) the choice of interpolation functions, and 
(3) the location of the grid points. The number of elements, and hence the 
number of grid points, is usually restricted simply by computer capability and 
by processing costs. Also, there have been many recent advances in improving 
the accuracy of the finite-element method by using higher order interpolation 
polynomials and shape functions (p-method) and by exhaustive analysis with 
large numbers of elements (h-method). Some of these advances are described in 
the references cited herein. However, optimal grid point location (r-method) 
is far less advanced. Practical procedures for the analyst still need to be 
developed and refined. 

One of the earliest attempts to develop a grid optimization procedure was 
that of Prager in 1975 (ref. 1). Prager's work provided a stimulus and a basis 
for later grid optimization research as recorded in references 2 to 19. Note- 
worthy among these efforts are the works of Shepard (refs. 11, 13, and 15), 
Masur (ref. 2), Turcke (refs. 8 and 9), Carroll (ref. 7), NcNiece (ref. 4), 
Carey (ref. 16), Diaz (ref. 12), Melosh (ref. 10), Durocher (ref. 17), and 
their col leagues. 

Prager examined a bar with a linearly varying cross section under tension. 
He showed that the grid producing the desired least potential energy is the one 
where the cross-section areas at the nodes form a geometric series. In this 
configuration, the strain energy is divided equally among the elements. 


Masur (ref. 2) observes that this latter result of equal element strain 
energies is not a general characteristic of optimal meshes but instead is a 
result of the simple geometry of Prager's problem. 

In this paper we present a finite-element grid improvement technique which 
is based on the minimization of the trace of the global stiffness matrix. We 
show that this method leads to identical results to those of Prager. It has 
the advantage of being simpler than traditional optimization procedures. 

The method presented herein provides a mesh improvement which is based on 
the geometry of the body. As such, it provides a significant improvement over 
uniform meshes, and it produces a good first iteration for accommodating spe- 
cial loading configurations. 

In the usual finite-element procedure, the governing equations are 
obtained by minimizing a functional ir by varying the dependent variables of 
the physical problem (ref. 20). For elastostatics this is equivalent to the 
principle of minimum potential energy (ref. 21). This leads to the familiar 
system of linear algebraic equations. Attempts to minimize ir with respect 
to the nodal coordinates, however, leads to a system of nonlinear equations. 
These equations are generally extremely difficult to solve even for the sim- 
plest cases. To avoid this difficulty we are proposing instead to examine the 
stiffness matrix to obtain information about an improvement in grid point loca- 
tion. Our motivation is the observation of the major role of the stiffness 
matrix in the value of the potential ir. Also, the entries of the stiffness 
matrix are dependent on the grid point coordinates. 


NOMENCLATURE 

A area 

A^ nodal area 

Aq, end area 

A 0 base area 

\ area of element 

c area ratio (see eq. (7)) 

E elastic modulus 
{f} global force array 

{f} transformed global force array 

f] entries of {f } 

E K ] global stiffness matrix 

EK ] diagonal form of CK3 
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entries of [K] 


L 

*k 

n 

P 


r ij 

S n 

T 

CT3 

u 

{u} 

{U} 

Ui 

X 

x k 

K 

ir 

5 

5k 


bar length 

length of k^* 1 element 
number of elements 
axial load 
length ratio, 0.^ /dj 

n 

sum of length ratios, E r., 

j=l 2 

trace of [K] 
transformation matrix 
axial displacement 
global displacement array 

transformed global displacement array 

/-'V 

entries of {u} 
axial coordinate 
nodal coordinate 
element stiffness matrix 
potential energy 
dimensionless axial coordinate 
dimensionless nodal coordinate 


ANALYSIS 

Our objective is to develop a practical and efficient procedure of grid 
enhancement tending towards optimization. Our thesis is that for many problems 
the minimization of the trace of the stiffness matrix leads to a minimization 
of the potential energy and, as a consequence, provides the optimal grid con- 
figuration. 

To see this, consider the governing matrix equation of finite-element 
analysis: 


EK]{u} = {f} (1) 

where [K] is the stiffness matrix, {u} the array of dependent variables, and 
{f } the force array. We can view [K] as an operator which maps {u} into 
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{f } . In this context, since [K] is symmetric 1 we can find an orthogonal 
transformation [T] which diagonalizes CK] ; that is, 

CK] = CT] t CK] CT] (2) 

where [K] is a diagonal matrix. Let [T] {u} and CT]{f} be {u} and {f } . 
Then the potential energy ir may be expressed as 

ir = \ {u} T [K]{u} - {f} T {u} = \ {u} T [K]{u} - {?} T {u} (3) 


In terms of the array components, ir becomes 


ir 


n 




i=l 



(4) 


/*>. /-s 

where the kj (i=l,...,n) are the diagonal entries of E K] . 


n ^ ^ 

Observe in equation (4) that the last term - E f.u. does not explicitly 

i =1 1 1 


involve the nodal coordinates. 


n _ ^ 
f i u i 


minimization of ir 


Therefore, - E 
i = l 

with respect to the nodal coordinates. 


does not affect the 


Also, since the u 


are positive and are independent variables in the minimization of ir, the mini- 
mization of ir with respect to the nodal coordinates occurs when the sum of 

the ki (the trace of [K]) is a minimum. Since the trace of a matrix is 


invariant under an orthogonal transformation, minimizing the trace of CK] is 
equivalent to minimizing the trace of CK]. 


In minimizing the trace, we will not adversely affect the diagonal domi- 
nance of CK] required to avoid ill-conditioning. The improved stiffness 
matrix we seek is the result of redistribution of the nodes and not of an arbi- 
trary mathematical operation. 

To illustrate the application of these concepts, consider the axially 
loaded tapered bar shown in figure 1. (This is the same problem examined by 
Prager (ref. 1) and Masur (ref. 2).) The objective is to determine a finite- 
element mesh which best predicts the axial displacement. Let the bar have 
length L and let it be divided into n elements with n + 1 nodes (numbered 

0 to n) as shown. Let the areas at the ends of the bar be A 0 and A&. Let 

1 be the nondimensional length parameter defined as 




x 

L 


( 5 ) 


^he analysis which follows is based on the symmetry of CK]. If [K] 
is not symmetric, a similar analysis could be developed using nonorthogonal 
transformations . 
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Then the area at any particular 5 along the bar is 


A = A 0 ( 1 - cO 


( 6 ) 


where c is 


c 



Hence, the area at the k th node is 


o < c < 1 


A|< = A 0 ( 1 - c£k> 


where Ck is ^ x k^> 


(7) 


( 8 ) 


Let the individual elements have a uniform_cross section. For example, 
let the k th element have cross-section area A^ and length &k as in fig- 
ure 2. (Note that the elements do not necessarily have the same length.) Then 
A|< and 8-k are 


A k = 


A + A 
k-1 k 




and 


> 


a k = x k " x k-l 


(9) 


- L( *k -<*-!> 

The element stiffness matrix for the k^h element is (ref. 20) 


A k E 

[ K] = f- 

K 


1 

-1 


-1 

1 


(10) 


where E is the elastic modulus. Then the trace T of the global stiffness 
matrix is 


T - 2E ) — = - 

1 ' fi.. L 


1 


' A i-1 + V 

Ji - «1-1. 


( 11 ) 


i = l 


i = l 


The improved node location occurs when the trace is minimized with respect 
to the nodal coordinates ^(k=l , . . . ,n-l ) .2 Hence, by setting the derivative 
of T with respect to 5k equal to zero, we obtain 


^Nodes numbers 0 and n are fixed at the ends of the bar and thus not 
considered for optimization. 
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E 

af; <\-i ♦ V 

3^ (A k * Vl* 

A k-1 + A k 

L 

*k - 5 k _i 

^k+1 ' ^k 

" - 5 k-i 


\ + A k+1 


’k+1 ^k 


Using equation (9) and simplifying we obtain 


A k + 1 = A k 


- 1 + A 


k-1 ft 


To simplify the analysis it is convenient to introduce the length ratio 
parameter r^j defined as ft^/dj. Then the ratio d|^ + i/Q.|^ may be written as 


Then L/2^ is 


k+1 ,1 
r.,, 


n n 

\ h \ hi.h 

r ki ' r ki 


^k L * k L 
j=l j=l 


where S is defined as E r.,. Using this notation, we can rewrite 
n j=l 31 

equation (13) as 


A = A [ JS+LJ 

k+1 k V r kl 


1 + A, 


k+1 , 1 
r. i 


c Vk + l 


Also, ^ may be written as 


ft l + a 2 + ‘ ‘ ’ + a k 


+ r 21 + r 31 + 
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Then, from equation (8), may be written as 


\ = A o 


1 - c 



n/J 


Specifically, A ] and A 2 are 


A 1 " A o 


and 


A 2 = A o 


1 - c 


n/J 


(1 + r 2 1 ) 




To obtain the element area ratios let k = 1 in equation <16). A , 


A, = A. r«,- 1 + A r!;, + cA r 


r 21 + 1 


2 ' "1 21 


o' 21 T V "V21 S 


Then, by using equations (19) to eliminate A 1 and A 2 , we obtain 


V 1 r 21 ) ■ c 


From the first of equations (19) we have 


or 


A 1 ‘ A o r 21 


□. _ r 

A 0 r 21 


Similarly, the second of equations (19) leads to 

2 


A 2 ” A o r 21 


or 


A, ~ r 21 






( 18 ) 

(19) 

is then 

(20) 

( 21 ) 

( 22 ) 

(23) 
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Next, let k = 2 in equation (16). By using the same procedure we obtain 


- r~„ = k— ( 1 - r.„ + r„, ) 


32 " S 


32 T '21 


But, in view of equation (21) this becomes 


or 


1 - r 32 - (1 - r 2 1 ) ( 1 r 32 + r 21 ) 


r 21 <r 32 " r 21 ) " 0 


•\ 




Thus, 


From equation (18) we have 


r 32 " r 21 


A- = 

3 O 


1 -elf 


- A r 3 
' A o r 21 


Therefore, we obtain 

A 3 A» A, 

AT = AT = A" = r 21 = r 32 

2 1 O 

Proceeding similarly for k = 3,4,.., we obtain 

A. Vi r r . r 

Vl ‘ \-2 A 0 ' 21 ' 32 V"-' 


Thus, we have the relations 


"\ 


r 31 ' r 2 1 r 32 “ r 21 


r 41 ~ r 43 r 32 r 21 = r 21 ’ ‘ ‘ ' ,r kl ~ r 21 


Hence, is the geometric series 


j 


1 - r 

2 k -1 21 

S k = 1 + r 2 i + r 2 i + ... + r 2] = (1 _ r^) 


(24) 


(25) 


(26) 


(27) 


(28) 


(29) 


(30) 


(31) 
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Then, from equations (18) and (21), we have 


\ - A o 




(1 - r^i )] = A_r 


.k 

21 


k 

o' 21 




and then 


A n “ A o r 21 " A 0, 




Then, from equation (7) we see that is 


( 32 ) 


r 2] = (1 - c) 


1/n 


(33) 


Finally, by substituting into equation (17) we have 


^k = L (1 + r 21 + r 31 + + r kl > _ (1 r 21 > 


1 + IT £ ^ + ^21 + • • • + ^ 


k-l 

21 


1 - r 


21 1 - (1 - c) 


k/n 


(34) 


This is the result obtained by Prager (ref. 1) in his analysis of the same 
problem. 


DISCUSSION 

First, observe that in equation (34) for a uniform thickness beam c = 0 
and thus ^ is undetermined. This means that for a uniform thickness beam 
the nodal positions are arbitrary; that is, all meshes are equally optimal for 
a uniform thickness beam. 

Next, consider again the element stiffness matrix of equation (10). From 
equations (8) and (34) the scalar multiplier is 


E( \-i - v v r * r 2i) 

l k = 2L<5 k - 5 k .,> ‘21- U - 


(35) 


Since this is a constant (independent of k) the element stiffness matrix is 
the same for each element. This means that each element has the same strain 
energy. Masur (ref. 2) has suggested that this result is due to the simple 
geometry of the problem. 

Even with this simple geometry, however, the analysis needed to determine 
the optimal nodal positions has been extremely detailed. With more complex 
geometries the analysis will become intractable. Alternatively, a more conven- 
ient method of improving the nodal positions is to examine the trace of the 
stiffness matrix and to adjust the nodal positions to minimize the trace. The 
criteria for minimizing the trace of the stiffness matrix is a comparatively 
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simple procedure - readily amenable to the development of computer algorithms 
for optimal nodal locations. 

In the preceeding discussion, we have ignored any consideration of the 
effects of concentrated loads or boundary/edge effects. It is common engineer- 
ing practice to use a fine (dense) grid in highly loaded regions and to place 
a grid point specifically at the point of application of a concentrated load. 
The trace minimization suggested here is intended to aid in the discretization 
process at locations removed from concentrated loads. Our justification is 
based on Saint-Venant 1 s principle that localized effects disappear at short 
distances . 

A principal question with the minimum trace method is: What is the range 

of applicability? The authors have found the method to lead to significant 
decreases in potential energy from that of a uniform mesh for structural and 
heat-transfer problems. The range of applicability is currently being 
explored. Numerical algorithms using this procedure are being developed. 
Finally, the influence of values of second and higher invariants of the stiff- 
ness matrix needs to be explored. 


NUMERICAL EXAMPLE 


To illustrate the value of optimizing the mesh, consider an axially loaded 
bar which tapers to 1/3 the base area as in figure 3. Specifically, let P, 

A 0 , c, E, and L have the following values: 


P = 20 N 

A = 0.0015 m 2 
o 



E = 2.07X10 11 N/m 2 
L = 4 m 




r ( 36 ) 




The objective is to find the axial displacement. 


From elementary mechanics the axial displacement u at any location x 
i s 


“ - - s jjh - f) <37> 

To compare the displacement results of finite-element models with equa- 
tion (37), four models of the bar, each having four elements, were examined. 

One of the models had a uniform nodal distribution. Another had the "optimal" 
mesh as developed in equation (34). The remaining two models had arbitarily 
selected nodal distributions. The nodal displacements were evaluated using the 
four models and compared with the displacement calculated by equation (37). 
Table I shows the results. Table II presents an error analysis and also an 
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l- 2 -norm of the error. As expected, the optimum mesh produces the least L 2 
error. 


CONCLUSIONS 

A new method of finite-element grid improvement based on the principle of 
minimization of the trace of the stiffness matrix was developed. This proce- 
dure is equivalent to minimizing the potential energy of the model by dividing 
the strain energy equally among the elements. The following conclusions are 
made: 


1. The analysis and the numerical results demonstrate the potential use- 
fulness of the trace minimization mesh improvement method. 

2. Minimization of the trace of the stiffness matrix is a relatively 
simple mesh optimization procedure. It is readily adaptable to algorithm 
development. 

3. Trace minimization can be used in combination with other grid optimiza- 
tion techniques. Indeed, it can provide a starting mesh for iteration tech- 
niques, useful for specialized loading. 

4. The principal benefit of the trace minimization method is accuracy of 
analysis as opposed to efficiency of analysis. 
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TABLE I. - COMPARISON OF AXIAL DISPLACEMENTS FOR THE BAR OF FIGURE 3 CALCULATED 

USING VARIOUS MODELS 


Axial 
location, 
x , 
m 

Exact displacement 
(eq. (37)), 

10 -9 m 

0.0 

0.0 

.5 

33.6276 

1 .0 

70.46244 

1.4409860 

106.1461 

2.0 

156.702 

2.5 

208.3078 

2.5358986 

212.2923 

3.0 

267.883 

3.3678522 

318.4384 

4.0 

424.5845 


Displacements computed using various models, 


Uni form 
mesh 


0.0 

70.2679 

156.1509 

266.5719 

421.1612 


"Optimum" 

mesh 

(Prager) 


0.0 

105.482 

210.9648 

316.4529 

421.940 


m 


Mesh 3 


0.0 

33.60639 

70.41338 


206.8158 


417.6195 


Mesh 4 


0.0 

33.60639 


15.56506 

207.1804 


417.9814 


Axi al 
locati on , 
x, 
m 


0.0 

.5 

1.0 

1.441 

2.0 

2.5 

2.536 

3.0 
3.368 

4.0 

L ^-Norm 


TABLE II. - ERROR ANALYSIS 


Error for various meshes, 
1(T 9 m 


Uni form 
mesh 


0.0 

.1945411 

.5506105 

1.311108 

3.423228 


"Optimum" 

mesh 

(Prager) 


0.0 

.664118 

1.32769 

1.985439 

2.64433 


Mesh 3 


0.0 

.02121 

.04906 


Mesh 4 


0.0 

.02121 


1.0514 

1.1274 


6.6031 


3.7119418 3.6246009 7.1232118 6.780697 


k TH ELEMENT 


FIGURE 1. - LINEAR TAPERED BAR 
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